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ABSTRACT.  This  paper  examines  the  scattering  effect  of  knit  lines  and  voids  in 
SOFI  through  simulations  of  THz  interrogation  at  normal  and  non-normal  angles  of 
incidence  and  using  focused  and  non-focused  single-cycle  plane  waves.  We  model 
the  electromagnetic  field  using  the  TE  mode  of  the  2D  Maxwell’s  equations  reduced 
to  a  wave  equation,  which  are  then  solved  with  a  finite-element  time-domain  method. 
The  knit  lines  are  modeled  by  changing  the  dielectric  constant. 
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INTRODUCTION 

Pulsed  THz  frequency  waves  have  been  shown  to  be  particularly  useful  for 
the  detection  of  voids  in  the  Sprayed  on  Foam  Insulation  (SOFI)  used  on  the 
Space  Shuttle’s  Thermal  Protection  System  (TPS)  [1].  However,  the  mod¬ 
eling  of  the  propagation  of  a  THz  pulse  inside  of  a  material  which  exhibits 
heterogeneous  micro-structures  of  sizes  that  are  on  the  order  of  the  wave 
length  of  the  interrogating  field  is  not  straight-forward.  Additionally,  data 
on  the  dielectric  properties  of  low  density  foam  in  the  THz  regime  is  rather 
sparse.  Initial  efforts  to  remedy  this  deficiency  can  be  found  in  [2]. 

As  discussed  in  [3],  usual  approaches  to  THz  interrogation  of  SOFI  gen¬ 
erally  employ  signal  processing  techniques  (for  example,  see  [4]),  which  do 
not  take  advantage  of  much  of  the  information  contained  in  the  reflected 
data  signal.  A  physics  based  model  may  be  able  to  more  accurately  describe 
defects  [5]. 

Previous  efforts  ([5,  6])  investigated  the  application  of  polarization  mech¬ 
anisms  to  account  for  the  attenuation  of  wave  propagation  in  foam,  but 
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matching  the  simulations  to  experimental  data  has  yet  to  be  completely 
successful.  While  these  models  were  expressed  in  only  one  spatial  dimen¬ 
sion,  the  2D  formulation  in  [3]  allowed  for  non-normally  incident  angles 
and  curved  interfaces.  Further,  by  allowing  the  dielectric  constant  to  be 
piecewise-constant,  layers  of  differing  densities  can  be  modeled  explicitly. 

SOFI  is  applied  in  layers,  producing  interfaces  between  layers  that  have 
increased  density  (called  “knit  lines”)  which  scatter  the  interrogating  wave¬ 
form.  As  the  knit  lines  are  generally  on  the  order  of  .5  mm  thick,  they  are 
frequently  ignored  in  simulations.  But,  as  there  may  be  as  many  as  20  knit 
lines  in  a  typical  20cm  block  of  foam,  the  aggregate  effect  on  the  interrogat¬ 
ing  signal  can  be  significant.  This  paper  describes  our  efforts  to  understand 
and  quantify  the  scattering  mechanisms,  as  well  as  to  minimize  their  effects 
by  the  use  of  focusing  and/or  altering  the  angle  of  incidence.  We  consider 
the  scattering  of  a  THz  plane  wave  off  of  the  knit  lines,  modeled  as  layers 
of  increased  density,  and  voids,  modeled  as  pockets  of  no  density,  inside  a 
block  of  low  density  polyurethane  foam. 

An  outline  of  the  paper  is  as  follows.  First  we  present  the  particular  form 
of  Maxwell’s  equations  we  use  to  model  the  2D  electromagnetic  interro¬ 
gation  problem  for  foam.  Then  an  experimentally  based  approach  for  esti¬ 
mating  the  material  parameters  for  use  in  the  system  is  described.  We  next 
briefly  report  on  our  numerical  techniques  for  solving  the  resulting  system. 
Finally  we  present  simulations  for  cases  in  which  the  angle  of  incidence  is 
altered,  and  the  plane  wave  is  focused  to  a  point. 

A  SIMPLIFIED  MATHEMATICAL  MODEL 

For  the  computational  domain  we  take  a  rectangular  region  0  <  x  <  0.1  m 
and  0  <  y  <  0.2  m.  Figure  1  shows  a  schematic  region  containing  a  material 
with  knit  lines  (represented  by  dashed  lines)  1mm  from  each  other,  each 
parallel  to  an  approaching  plane  wave,  and  perpendicular  to  the  direction 
of  propagation.  The  regions  between  knit  lines  consist  of  a  low  density 
material,  and  the  elliptical  pocket  (void)  is  modeled  as  having  zero  density. 
The  far  right  boundary  (y  =  0.2)  is  assumed  to  be  metallic,  and  therefore 
supra-conducting,  thus  simulating  the  aluminum  backing  of  the  SOFI  on  the 
shuttle  external  tank.  A  vacuum  is  present  to  the  left  of  the  material  in  which 
the  interrogating  field  will  be  generated. 
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FIGURE  1:  Schematic  of  a  plane  wave  approaching  a  domain  containing  knit  lines  and  an 
elliptical  void.  Dashed  lines  represent  knit  lines,  dot-dash  is  foam/air  interface.  Elliptical 
pocket  (5  mm  wide)  between  knit  lines  is  a  void  modeled  by  n  =  1.  “+”  marks  the  signal 
receiver.  Back  wall  is  perfect  conductor. 

We  combine  the  TE  mode  of  the  two  dimensional  Maxwell’s  equations 
into  one  equation 

+ v-  (1) 

where  e{x)  and  fi(x)  are  the  spatially  dependent  dielectric  permittivity  and 
permeability,  respectively.  The  corresponding  speed  of  propagation  is 

1 

e{x)ii{x)  ’ 

where  cq  is  the  speed  in  a  vacuum. 

For  our  source  current,  Js,  we  want  to  simulate  a  windowed  pulse,  in  this 
case  a  pulse  that  is  allowed  to  oscillate  for  one  half  of  one  period  and  then 
is  truncated.  Although  generators  produce  a  curved,  sometimes  spherical 
wave,  we  assume  the  target  is  sufficiently  far  from  the  generator  (approxi¬ 
mately  6”)  so  that  the  wave  is  essentially  planar  when  it  reaches  our  domain 
of  interest.  Thus,  we  want  the  pulse  to  originate  at  x  =  0,  the  beginning  of 


3 


our  computation  domain,  which  we  model  as  a  delta  distribution  centered  at 
x  =  0.  In  order  to  have  a  smooth  source  we  use  a  function  of  the  form 

Js(t,x )  =  S(x)e~^t~t°^t°^  ,  (2) 

where  to  =  f//4  when  tf  is  the  period  of  the  interrogating  pulse.  For  exam¬ 
ple,  if  the  frequency  is  /  =  .2THz,  then  tf  =  1/ /  =  1  x  10_11s.  A  typical 
value  for  the  exponent  is  b  =  4. 


Boundary/Initial  Conditions 


To  model  a  metallic  backing  behind  the  foam,  we  use  reflecting  (Dirichlet) 
boundary  conditions 

[E]x=o.\  =  0. 

In  order  to  have  a  finite  computational  domain,  we  impose  first  order  ab¬ 
sorbing  boundary  conditions  at  x  =  0,  which  are  modeled  as 


dE  ,,dE 


=  0. 

x=0 


With  these  boundary  conditions,  ideally  a  normally  incident  signal  passes 
out  of  the  computational  domain,  and  does  not  return,  i.e.,  we  force  it  to 
be  absorbed  by  the  boundary.  Note  that  for  signals  that  are  incident  at  an 
angle,  some  reflection  occurs.  Lastly,  to  allow  for  propagation  along  the  top 
and  bottom  boundaries  (y  =  0  and  y  =  0.2),  we  use  insulating  boundary 
conditions 


'dE' 

=  o, 

'dE' 

_dy  _ 

y= o 

.  dy  _ 

Also  we  assume  zero  initial  conditions,  i.e., 


^(O,  x)  =  0,  ^(O,  x)  =  0. 


Modeling  Knit  Lines 

The  speed  of  propagation  in  the  domain  is  given  by 

1 

e(x)y  o’ 
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where  Co  is  the  speed  in  a  vacuum  and  n  is  the  index  of  refraction.  We 
model  knit  lines  by  changing  the  index  of  refraction,  thus  effectively  the 
speed,  in  that  region.  In  order  to  model  the  speed  of  wave  propagation 
in  the  knit  lines  versus  the  material  surrounding  them,  we  need  to  be  able 
to  distinguish  between  the  respective  indices  of  refraction,  i.e.,  n\  in  the 
low  density  region,  and  rio  in  the  higher  density  knit  line.  Further,  we  can 
currently  only  measure  the  effective  index  of  refraction  of  the  composite 
material,  ne,  which  is  done  by  computing  the  “time  of  flight”  in  experiments. 
Thus  we  need  to  relate  these  three  indices  to  each  other  in  order  to  have 
accurate  estimates  for  the  propagation  speed  in  each  region,  to  use  in  the 
simulations. 

Experiments  have  suggested,  via  time-of-flight  measurements,  a  value 
for  the  effective  index  of  refraction  as  ne  =  1.032254=0.001.  To  estimate  the 
index  of  refraction  for  foam  in  the  absence  of  knit  lines,  blocks  of  5 cm  x  5cm 
blocks  of  foam  with  varying  numbers  of  knit  lines  layers  (between  3  and 
7,  approximated  to  the  nearest  .25)  were  produced  and  interrogated.  By 
calculating  the  velocity  of  the  pulse  in  these  blocks  we  may  extrapolate  the 
data  to  the  case  of  zero  knit  line  layers.  Figure  2  displays  the  data  collected 
in  the  experiment. 

A  line  of  best  fit,  computed  using  linear  regression,  is  plotted  in  Figure  2 
with  error  bars  of  two  standard  deviations.  The  range  of  values  correspond¬ 
ing  to  zero  knit  lines  is  the  estimate  for  the  velocity  of  the  pulse  in  a  region  of 
foam  without  knit  lines.  The  mean  value  in  this  region  is  2.94638  x  108m/s 
resulting  in  an  index  of  refraction  for  the  low  density  region  ofni  =  1.0172. 
Using  ne  =  (1  —  v)n\  +  /vn2,  where  v  is  the  volume  fraction  of  the  knit 
line  region  in  the  foam  block.  The  volume  fraction  v  can  be  estimated  by 
noting  the  thickness  of  each  knit  line  divided  by  the  period  in  which  the  knit 
lines  occur.  Assuming  0.5mm  knit  lines  in  each  1cm  of  foam  corresponds 
to  v  =  .05.  This  gives  the  index  in  the  knit  line  as  n2  =  1.1869. 

In  [3],  the  effective  index  of  refraction  ne  was  approximated  via  the  clas¬ 
sical  Claussius-Mossotti  equation  (for  instance,  see  [7]  and  [8]),  by  assum¬ 
ing  the  total  polarizability  is  the  sum  of  the  two  polarizabilities  in  each  part 
of  the  material.  The  index  of  refraction  in  each  part  of  the  foam  was  esti- 
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FIGURE  2:  Experimental  data  for  the  velocity  of  a  THz  pulse  versus  the  number  of  knit 
lines,  with  the  line  of  best  fit  from  applying  linear  regression.  The  error  bars  represent  two 
standard  deviations;  the  “X”  denotes  the  Claussius-Mossotti  estimate. 

mated  to  be 


m  =  1.01398,  (3) 

n2  =  1.03507.  (4) 

This  value  is  plotted  as  an  “X”  in  Figure  2,  and  is  within  the  range  of  the 
error  bars  of  the  linear  regression  estimate. 

Although  the  experimental  result  agrees  well  with  the  theoretical  esti¬ 
mate,  it  is  clear  from  the  scattered  nature  of  the  data  points  in  Figure  2  that 
the  correlation  between  the  number  of  knit  lines  and  the  velocity  is  not  per¬ 
fectly  linear.  The  plot  in  Figure  3  shows  that  this  is  mainly  due  to  the  rather 
imprecise  relationship  between  the  density  and  the  number  of  knit  lines. 
The  symbol  in  the  plot  denotes  the  frequency  (number  of  occurrences)  of 
data  point.  The  solid  line  again  represents  a  line  of  best  fit.  Only  one  mea¬ 
surement  of  the  density  is  beyond  two  standard  deviations  from  the  mean, 
but  many  are  outside  of  one  standard  deviation.  Better  techniques  in  esti¬ 
mating  the  number  of  knit  lines  in  a  block  of  foam  may  result  in  density 
data  more  closely  linearly  related  to  the  number  of  knit  lines,  and  therefore 
a  more  precise  estimate  of  the  density  for  foam  in  the  limit  of  zero  knit  lines. 
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FIGURE  3:  Experimental  data  for  the  density  of  a  5cm  x  5cm  block  versus  the  number  of 
knit  lines,  with  the  line  of  best  fit  from  applying  linear  regression.  The  error  bars  represent 
two  standard  deviations;  the  symbol  denotes  the  frequency  (occurrence)  of  each  data  point 
present. 


NUMERICAL  METHOD 

We  employ  a  (second  order)  Finite  Element  method  using  standard  linear 
two  dimensional  (Q 1)  basis  elements  to  spatially  discretize  the  model  de¬ 
scribed  by  (1).  This  results  in  nine-banded  mass  and  stiffness  matrices,  M 
and  S.  We  also  have  a  contribution  from  absorbing  boundaries  which  we 
denote  with  B.  Thus  our  semi-discrete  system  for  the  vector  of  electric  field 
values  e  is 

Me  +  Be  +  Se  =  f. 

Here  we  have  absorbed  coefficients  \  and  -  into  the  definitions  of  M  and 

Cz  C 

B,  respectively. 

For  the  time  derivatives  we  use  second  order  discretizations  (centered 
differences)  for  both  the  first  and  second  derivatives.  After  collecting  all 
terms  involving  the  updated  time  step  into  the  left  hand  side  of  the  equation 
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we  have  the  following  linear  system 

-A6n-i-i  =  6,  (5) 

where  A  contains  multiples  of  M  and  B,  and  b  depends  on  en  and  en-\,  as 
well  as  S  and  /. 

As  discussed  in  [3],  a  mass-lumping  approach,  using  a  quadrature  rule 
applied  to  the  basis  functions,  is  most  efficient  for  solving  the  linear  sys¬ 
tem  (5)  on  each  time  step  considering  the  length  of  time  span  required  for 
propagation  problems  of  this  type.  LU  factorization  is  limited  by  memory 
constraints  for  large  problems,  and  iterative  methods  are  restricted  by  com¬ 
putational  time. 

SIMULATIONS 

We  compute  numerical  simulations  of  an  electromagnetic  wave  propagating 
through  a  material  described  by  its  index  of  refraction,  which  determines  the 
speed  of  propagation.  We  consider  the  presence  of  a  void  similar  to  what  is 
seen  in  SOFI  when  a  layer  does  not  completely  fill  a  recess  formed  in  the 
previous  layer.  As  the  foam  cures,  an  elliptical  pocket  of  air  is  formed.  The 
void  is  modeled  by  taking  its  index  of  refraction  to  be  that  of  free  space, 
i.e.,  no  =  1.  The  knit  lines  are  modeled  with  index  of  refraction,  77.2,  and 
the  surrounding  low  density  regions  are  described  by  n\.  Throughout  the 
following,  the  knit  lines  are  described  using  the  values  estimated  above  by 
the  use  of  experimental  data. 

Plane  Wave 

Figure  4  depicts  snapshots  in  time  of  the  propagation  of  a  plane  wave  inci¬ 
dent  on  a  void  in  the  material.  The  cross  denotes  the  location  of  the  receiver 
to  collect  data.  The  reflection  from  the  void  is  clearly  seen  in  the  third 
frame.  This  reflection  expands  out  to  form  an  oblong  elliptical  wave  which 
eventually  returns  to  the  receiver  where  the  signal  is  recorded. 

As  the  actual  transceiver  in  the  experimental  setup  is  devised  of  a  trans¬ 
mitter  and  receiver  which  do  not  occupy  the  same  space  (see  Figure  5),  the 
direction  of  pulse  propagation  is  at  some  non-normal  angle  to  the  plane  of 
the  foam  and  the  metallic  backing.  The  simulations  of  the  scenario  where 
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FIGURE  4:  Surface  plots  of  solutions  for  the  case  where  the  plane  wave  pulse  is  normally 
incident  on  the  material. 


FIGURE  5:  Picometrix  T-Ray  setup  used  in  the  electromagnetic  interrogation  of  foam. 
Note  the  non-normal  incidence  and  ability  to  focus. 
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the  plane  wave  pulse  is  incident  on  the  material  at  a  non-normal  angle  are 
displayed  in  Figure  6.  While  the  receiver  still  located  at  x  =  0.03,  it  has 
been  raised  in  the  y  direction  to  collect  the  center  of  plane  wave  reflection. 


FIGURE  6:  Surface  plots  of  solutions  for  the  case  where  the  plane  wave  pulse  is  non- 
normally  incident  on  the  material. 


The  plots  in  Figure  7  display  the  simulated  data  collected  at  the  receiver 
for  normal  and  non-normal  angles  of  incidence.  In  the  plot  the  relative  mag¬ 
nitude  of  the  reflection  versus  the  interrogating  signal  is  apparent.  The  inset 
displays  a  plot  of  the  magnification  of  the  reflection  from  the  void.  There  is 
a  not  a  significant  difference  in  the  two  reflections,  only  a  slight  timing  dis¬ 
crepancy.  Note  that  nearly  all  of  original  signal  returns  even  with  an  oblique 
angle  of  incidence. 

Focused  Wave 

For  a  focused  source  we  model  the  interrogating  signal  using  scattered  field 
formulation  of  point  source  reflected  from  elliptical  mirror.  Note  that  the  top 
and  bottom  boundary  conditions  are  now  absorbing  as  we  no  longer  have  a 
plane  wave  propagating  along  these  boundaries. 
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Signal  received  at  x=0.03  Signal  received  at  x=0.03 


FIGURE  7:  The  left  (right)  plot  displays  the  signal  received  at  x  =  0  on  center  line  for  the 
normally  (non-normally)  incident  plane  wave  case  (the  inset  plot  is  a  magnification  of  the 
reflection  from  the  void). 


Figure  8  displays  snapshots  in  time  of  the  simulated  propagation  of  a 
focused  wave  incident  on  a  void  in  the  material.  The  reflection  from  the 
void  is  clearly  seen  in  the  second  frame.  This  reflection  expands  out  to  form 
an  oblong  elliptical  wave  which  eventually  returns  to  the  antenna  where  the 
signal  is  recorded  with  a  receiver. 

The  simulations  of  the  scenario  where  the  focused  wave  pulse  is  incident 
on  the  material  at  a  non-normal  angle  are  displayed  in  Figure  9.  The  source 
is  modeled  using  scattered  field  formulation  of  point  source  reflected  from 
elliptical  mirror.  The  receiver  is  located  at  x  =  0.03,  but  raised  to  collect 
center  of  the  focused  wave  reflection. 

The  plots  in  Figure  10  give  the  simulated  data  collected  at  the  receiver 
for  the  focused  wave  incident  at  the  two  angles  of  incidence,  respectively. 
The  inset  again  displays  a  magnified  plot  of  the  reflection  from  the  void. 

There  is  a  distinct  difference  in  the  structure  of  the  two  reflections.  In 
particular,  the  front  part  of  the  void  reflection  is  clearly  more  pronounced  in 
the  normally  incident  case.  This  suggests  that  simulations  which  do  not  take 
the  proper  angle  of  incidence  into  account  may  grossly  over-exaggerate  the 
magnitude  of  the  reflection  from  the  void.  Note  that  the  amplitude  of  the 
void  reflections  in  both  cases  is  larger  than  that  of  the  plane  wave  scenario. 
Further,  the  difference  between  the  reflection  from  the  beginning  of  the  void 
is  more  distinct  from  the  reflection  off  the  back  of  the  void  for  the  focused 
wave.  As  the  energy  from  the  interrogating  wave  is  concentrated  on  the 
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FIGURE  8:  Surface  plots  of  solutions  for  the  case  where  the  focused  wave  pulse  is  nor¬ 
mally  incident  on  the  material. 


void  area,  and  thus  reflected,  more  information  is  returned  to  be  collected 
by  the  receiver.  Conversely,  although  reflections  off  the  void  are  larger,  the 
total  energy  that  returns  is  less  than  the  plane  wave  simulation.  While  these 
facts  about  focusing  are  not  new  results,  it  is  important  that  any  model  of  a 
focused  pulse  validate  these  basic  principles. 

CONCLUSIONS 

We  have  developed  a  framework  which  accounts  for  the  presence  of  knit 
lines  in  modeling  the  electromagnetic  propagation  of  interrogating  pulses 
in  SOFI.  We  were  able  to  compute,  using  results  from  a  laboratory  experi¬ 
ment,  estimates  for  the  index  of  refraction  in  both  the  knit  lines  and  the  sur¬ 
rounding  low-density  foam  based  on  measured  values  from  time-of-flight 
experiments  and  observable  material  properties.  The  values  were  shown  to 
be  consistent  with  classical  physics-based  theory  assuming  constant  polar- 
izibility. 

In  the  effort  to  detect  flaws  in  low-density  materials  such  as  foam,  highly 
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FIGURE  9:  Surface  plots  of  solutions  for  the  case  where  the  focused  wave  pulse  is  non- 
normally  incident  on  the  material. 


Signal  received  at  x=0.03  Signal  received  at  x=0.03 


FIGURE  10:  The  left  (right)  plot  displays  the  signal  received  at  x  =  0  on  center  line  for 
the  normally  (non-normally)  incident  focused  case.  The  inset  plot  is  a  magnification  of  the 
reflection  from  the  void. 


accurate  models  must  exist  to  give  simulations  the  precision  necessary  to 
distinguish  small  amplitude  reflections  from  noise,  including  from  model 
error.  While  the  current  formulation  may  be  used  as  a  forward  solution  in 
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an  inverse  problem  context,  it  is  likely  that  the  highest  value  will  lie  in  its 
ability  to  generate  synthetic  data  with  which  to  test  faster  signal  processing 
approaches  to  damage  detection.  Thus  it  can  be  used  either  to  explore  which 
shapes  of  voids  are  the  hardest  to  detect,  or  to  generate  data  for  scenarios 
that  are  difficult  or  expensive  to  manufacture. 

Our  efforts  here  provide  an  approach  to  enhance  the  accuracy  of  a  model 
by  making  it  more  representative  of  the  material  in  question  and  the  sig¬ 
nal  used  to  interrogate  that  material.  The  results  themselves  suggest  that 
accurate  modeling  of  the  angle  of  incidence  is  important  to  quantify  the 
amplitude  of  the  reflected  signal,  while  modeling  the  focusing  captures  the 
clear  distinction  between  the  reflections  from  the  front  versus  the  back  of 
the  void. 

Sufficient  experimental  results  for  comparison  are  difficult  to  obtain  since 
the  amplitude  of  reflections  off  of  low  density  materials  is  very  low  to  im¬ 
measurable  using  currently  available  power  sources  at  the  THz  frequency. 
Therefore,  the  results  of  this  paper  should  be  considered  as  motivation  for 
the  development  of  high  power  THz  devices.  For  now,  experimental  infor¬ 
mation  is  collected  from  the  aluminum  reflection.  More  work  needs  to  be 
done  to  match  simulations  to  this  data,  including  adding  attenuation,  possi¬ 
bly  via  coupling  Maxwell’s  equations  to  a  time  domain  model  of  a  scattering 
mechanism. 
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